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1.  Introduction 


Assessment,  delineation,  and  treatment  of  traumatic  brain  injury  (TBI)  are  critically  needed 
across  civilian  and  military  populations.  The  national  cost  of  TBI  is  estimated  to  be  $60  billion 
annually,  with  3.2  million  Americans  living  with  disabilities  due  to  TBI  ( 1 ).  Blast-related  TBI  is 
of  particular  interest  due  to  its  prevalence  in  recent  military  conflicts  (2).  Blast  exposure  results 
in  a  complex  loading  that  presents  many  possible  mechanisms  of  injury.  Diffuse  axonal  injury  is 
one  particular  injury  mechanism  that  has  been  cited  as  a  predominant  signature  injury  of  TBI 
neural  damage  (2,  3).  This  injury  is  characterized  by  widespread  structural  lesions  in  white- 
matter  fiber  tracts,  the  axons  of  neurons.  These  tracts  connect  brain  regions  into  a  structural 
network  and  allow  neurons  to  communicate  with  one  another  (4,  5).  Degraded  structural 
connectivity  has  been  linked  to  disease  states  ( 6 ,  7),  and  it  may  underlie  the  cognitive  deficits 
characteristic  of  mild,  moderate,  and  severe  cases  of  TBI  (S).  This  report  presents  a 
multidisciplinary  modeling  effort  that  aims  to  improve  our  understanding  of  how  mechanical 
loading  to  the  head  is  related  to  changes  in  the  structural  network  of  the  brain. 

This  modeling  effort  sits  within  a  larger  research  program  at  the  U.S.  Army  Research  Laboratory 
(ARL)  to  understand  brain  structure-function  couplings.  The  overall  goal  of  this  effort  is  to  better 
understand  how  variability  in  structural  connectivity  due  to  traumatic  insult  or  natural  variability 
between  healthy  individuals  relates  to  differences  in  the  brain’s  functional  connectivity  (9, 10) 
and,  ultimately,  to  individual  differences  in  human  behavior.  The  modeling  effort  described  here 
focuses  on  understanding  the  effects  of  traumatic  brain  insult.  The  use  of  a  finite  element  (FE) 
simulation  to  degrade  the  structural  brain  network  was  first  described  in  Kraft  et  al.  (11).  This 
report  aims  to  make  improvements  to  the  structural  network  weighting  and  degradation  methods. 
We  examine  if  the  graph  metrics  properties  used  to  characterize  the  damaged  structural  network 
are  scale-invariant  to  show  damage  consistently  across  different-sized  networks.  These 
improvements  are  needed  to  enable  future  research  efforts  that  will  couple  the  simulated 
structurally  damaged  network  with  a  more  coarse-grained  electrophysiological  model  developed 
by  David  and  Friston  (12)  and  implemented  at  ARL  (9, 13, 14).  The  electrophysiological  model 
simulates  functional  electrical  oscillations  for  a  simulated  brain  region  (node)  in  a  defined 
structural  network.  The  neural  mass  at  each  node  represents  a  small  patch  of  the  brain  containing 
several  thousand  interconnected  neurons.  Our  research  effort  provides  a  unique  capability  to 
examine  how  functional  ramifications  (e.g.,  abnormal  electrical  oscillations  between  brain 
regions)  are  tied  to  parameters  of  blast-induced  traumatic  brain  insults  by  using  the  damaged 
structural  network  from  an  FE  simulation  as  input  in  the  neurophysiological  modeling.  The 
long-term  vision  is  to  link  how  changes  in  these  electrical  oscillations  interfere  with  functional 
connectivity  patterns  that  are  characteristic  of  successful  performance  on  behavioral  tasks 
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(15-20).  Thus,  the  integration  of  two  modeling  efforts  with  ongoing  experimental  efforts  provide 
a  technical  framework  to  investigate  the  well-known  but  poorly  understood  cognitive 
impairments  that  often  result  from  traumatic  brain  insults  (9, 10). 

The  modeling  research  described  in  this  report  makes  use  of  concepts  from  structural  FE 
analysis,  structural  neuroimaging  methodologies,  and  graph  theory  analysis  approaches  from 
network  science.  In  short,  these  tools  are  combined  to  improve  the  simulated  mechanical 
response  of  brain  tissue,  estimate  resulting  damage  to  the  structural  connectivity  in  the  brain,  and 
develop  graph-theoretic  analysis  approaches  to  quantify  the  resulting  changes  in  the  structural 
network.  The  methods  described  incorporate  biological  thresholds  for  brain  tissue  damage  that 
are  used  to  degrade  the  brain’s  structural  connections  and  compute  a  resulting  damaged 
structural  network. 

The  FE  method  enables  a  physics-based  computational  simulation  of  a  mechanical  loading.  The 
strength  of  this  method  is  that  it  provides  the  full  three-dimensional  mechanical  response  of  a 
material  in  time,  thus  providing  insight  into  stress  and  strain  that  might  not  be  possible  to 
measure  in  vivo.  Finite  element  simulations  have  been  performed  for  both  impact  (21-29)  and 
blast  loading  (30-32)  to  capture  the  mechanical  response  of  the  human  head  and  brain.  However, 
these  investigations  do  not  directly  incorporate  their  results  into  plausible  models  of  injury  on  a 
network  level.  The  approach  to  link  physics-based  mechanical  parameters  and  biological  injury 
is  an  open  research  question.  Our  innovation  is  the  combination  of  the  FE  method  with 
neuroimaging  to  incorporate  biological  information  into  the  model  and  improve  mechanical 
response  estimates  and  biological  fidelity  of  cellular  injury  estimate. 

A  structural  neuroimaging  method,  known  as  diffusion- weighted  imaging,  allows  in  vivo 
imaging  of  white-matter  fiber  tracts  that  are  widely  implicated  in  long-distance  communication 
in  the  brain.  This  technique  uses  a  magnetic  resonance  imaging  (MRI)  scanner  to  image  the 
diffusion  of  water  in  the  brain.  This  directional  movement  reveals  the  local  brain  structure  since 
water  diffuses  in  the  same  direction  as  the  local  fiber  tract.  Postprocessing,  reconstruction 
tractography  algorithms  then  estimate  an  individual’s  whole-brain  fiber  tract  structure.  These 
fiber  tracts  constitute  the  structural  network  of  the  brain.  It  is  possible  to  then  characterize  this 
network  using  analytic  approaches  from  network  science. 

Graph  theory  analysis  methods  from  network  science  quantify  topological  properties  of 
networks,  and  they  have  been  successfully  adopted  for  characterization  of  brain  connectivity 
(33).  The  brain  is  divided  into  smaller  regions  of  interest  (ROIs)  to  provide  graph  nodes,  and  the 
white-matter  connections  between  the  ROIs  provide  the  graph  edges.  Standard  measures  of 
graph  properties  can  then  be  computed  to  characterize  the  structural  brain  network.  We  focus  on 
seven  key  measures  in  our  analysis:  degree,  shortest  path,  efficiency,  betweenness,  clustering 
coefficient,  modularity,  and  assortivity.  The  degree  measure  is  computed  for  each  node,  which  is 
the  sum  of  edge  weights  connected  to  a  node.  A  short  path  in  the  network  is  the  number  of  edges 
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traversed  or  the  sum  of  weights  along  the  most  direct  route  from  one  node  to  another.  Efficiency 
can  be  calculated  as  a  nodal  measure  referred  to  as  local  efficiency,  or  for  the  network  as  a  whole 
as  global  efficiency.  This  value  represents  how  well  connected  the  network  is  compared  with  a 
perfectly  connected  network  (34)  and  captures  the  networks  capacity  for  communication  along 
short  paths  (55).  Betweenness  is  the  number  of  times  a  node  is  part  of  the  shortest  path  between 
two  other  nodes  (34).  A  high  measure  of  betweenness  is  an  indication  of  a  hub  used  to  allow 
efficient  long-distance  communication  of  nodes  (34).  The  clustering  coefficient  is  a  measure  of 
how  well  a  particular  node  and  its  neighbors  are  connected  to  each  other  (34).  Modularity  is  a 
value  that  can  indicate  the  presence  of  a  community  structure  composed  of  localized  clusters  of 
nodes  (36).  The  assortativity  quantifies  how  similar  nodes  are  to  each  other  in  terms  of  their 
degree  (34). 

The  graph  theoretic  framework  has  been  applied  to  examine  effects  of  lesions  on  network  models 
of  the  brain.  Several  groups  have  investigated  the  effects  of  lesions  on  network  models  of  the 
brain  in  both  animal  models  (37,  38)  and  humans  (7,  39).  Young  et  al.  (37)  provide  an  early 
example  of  applying  damage  to  a  structure-function  coupled  network  based  on  the 
thalamocortical  system  of  a  cat.  They  use  anatomical  data  to  create  a  weighted  structural  network 
and  simulate  the  functional  network.  A  lesion  is  simulated  by  removing  a  node  from  the  network 
to  show  that  there  are  functional  effects  that  extend  beyond  the  lesion  site  (37).  Alstott  et  al.  (39) 
extend  this  concept  to  the  human  brain  by  using  diffusion-weighted  imaging  to  create  a  whole 
brain  structural  network.  They  show  that  specific  structural  network  measures,  such  as 
betweenness  and  efficiency,  can  characterize  the  network  and  quantify  the  effect  of  removing 
nodes.  These  research  efforts  establish  the  use  of  networks  to  describe  the  brain  as  a  whole  as 
well  as  with  simulated  lesions.  Our  work  will  complement  these  efforts  by  developing  injury 
models  that  degrade  the  network  to  represent  specific  loading  conditions. 

We  present  a  combination  of  the  FE  method  with  neuroimaging  and  network  modeling  to 
incorporate  biological  information  into  the  model  and  provide  the  capability  to  evaluate  the 
global  effect  of  specific  blast  loading  conditions.  The  work  presented  in  this  report  represents 
several  accomplishments  that  expand  on  previous  work  (11).  In  this  work  we  discuss  the  edge¬ 
weighting  method  and  how  it  is  modified  to  make  use  of  tractography  fiber  segments  rather  than 
whole  fibers.  This  prevents  the  complete  removal  of  edges  from  the  structural  network  seen  in 
previous  work.  Sensitivity  of  the  network  is  then  evaluated  by  comparing  two  simulations  that 
represent  a  frontal  blast  loading  and  a  blast  loading  from  the  side  to  show  the  ability  to 
differentiate  between  these  loading  directions.  Finally,  the  effect  of  scaling  the  number  of  ROIs 
in  the  structural  network  is  investigated.  An  evaluation  of  network  measures  is  performed  to 
determine  which  are  scale-invariant  between  networks  of  different  sizes  after  damage  is  applied. 
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2.  Methods 


We  briefly  summarize  some  of  our  methods  of  this  research  and  organize  them  as  follows. 
Section  2.1  describes  the  procedure  used  to  create  the  FE  mesh  from  neuroimaging  data.  Section 
2.2  provides  detail  of  the  FE  model,  including  the  loading  conditions,  material  properties,  and 
results  of  the  simulation.  Section  2.3  shows  how  FE  results  are  formatted  for  use  with  the 
structural  network  and  brain  atlas.  In  section  2.4,  we  introduce  our  new  fiber  segment  network 
degradation  method  and  discuss  how  it  is  different  from  Kraft  et  al.  (11). 

2.1  Neuroimaging  and  FE  Mesh  Creation 

Methods  to  create  the  FE  model  and  reconstruct  data  from  diffusion- weighted  imaging  are 
described  in  detail  in  Kraft  et  al.  (11,  40)  and  only  summarized  here.  In  short,  geometry  is 
created  from  MRI  data  by  using  the  software  Amira  (41)  to  segment  the  head  into  white  matter, 
gray  matter,  cerebrospinal  fluid,  skull,  skin,  and  other  soft  tissue  outside  of  the  skull.  In  addition 
to  the  head,  a  torso  from  the  Open  3D  Project  (42)  is  included  to  more  accurately  capture  blast 
propagation  to  the  brain.  This  geometry  is  used  to  generate  a  tetrahedral  mesh  for  FE  simulation. 
Tractography  is  reconstructed  from  a  diffusion-weighted  imaging  technique  called  diffusion 
tensor  imaging.  This  technique  provides  information  of  diffusion  as  a  tensor  value  that  can  be 
used  to  reconstruct  the  fiber  tractography  (43). 

2.2  FE  Simulation 

Here  we  briefly  summarize  some  of  the  key  features  of  the  model;  details  of  our  material 
parameters  and  constitutive  equations  are  described  in  detail  in  Kraft  et  al.  (11,  40). 
Tractography  from  diffusion  tensor  imaging  data  is  incorporated  into  the  three-dimensional  FE 
model  using  a  transverse  isotropic  material  model  specifically  developed  for  representing  white- 
matter  tissue  (44,  45)  where  each  FE  is  assigned  an  orientation  based  on  the  superimposed 
tractography  (40).  Material  constitutive  descriptions  and  properties  for  all  components  are 
described  in  detail  in  Kraft  et  al.  (11).  These  material  descriptions  include  properties  for  the 
skull,  cortex,  brain  stem,  cerebrospinal  fluid,  and  soft  tissue,  which  is  an  homogenized  mixture 
of  muscle  and  skin.  The  volume  of  the  body  is  also  modeled  as  an  homogenous  soft  tissue  to 
more  accurately  capture  pressure  loading  that  is  transmitted  to  the  head.  A  nodal-based 
tetrahedral  formulation  is  used  to  resolve  locking  issues  that  are  related  to  standard  tetrahedral 
elements  to  better  represent  soft  tissues  (46). 

The  transient  response  of  the  blast  loading  to  the  head  is  simulated  in  Sierra  (Sandia  National 
Laboratories),  an  explicit  FE  method  Lagrangian  code.  The  blast-loading  method  combines  the 
TNT  air  blast  work  of  Kingery  and  Bulmash  (47)  and  Randers-Pehrson  and  Bannister  (48),  with 
Sachs  scaling  to  match  the  ConWep  code  (49).  The  code  accounts  for  the  angle  of  incidence  by 
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transitioning  from  reflected  to  incident  pressure  by  the  relationship  described  in  Gullerud  et  al. 
(46).  For  the  comparison  examined  in  this  paper,  two  blast  conditions  are  simulated.  A  frontal 
loading  is  simulated  by  a  TNT  spherical  air-blast  equivalent  of  7.5  lb  located  0.4  m  above  the 
ground  and  a  2-m  standoff  distance  in  the  front.  A  side  loading  is  simulated  by  a  TNT  spherical 
air-blast  equivalent  of  7.5  lb  located  0.4  m  above  the  ground  and  a  2-m  standoff  distance  from 
the  right  of  the  body. 

The  results  of  the  two  loading  conditions  are  presented  in  figure  1  as  time  history  plots  of  the 
pressure,  axonal  strain,  and  strain-rate  for  the  frontal  loading  (panels  A-C)  and  side  loading 
(panels  D-F).  Pressure  is  positive  in  compression.  Axonal  strain  is  defined  as  the  strain  in  the 
direction  of  the  axon  based  on  the  tractography  within  an  element.  Strain  rate  is  an  effective 
measurement  based  on  the  strain-rate  tensor.  These  parameters  have  been  used  to  predict  cellular 
injury  within  neuronal  cellular  cultures  (50).  The  five  regions  of  interest  plotted  include  the 
frontal  lobe  (red),  temporal  lobe  (green),  occipital  lobe  (blue),  parietal  lobe  (orange),  and  corpus 
callosum  (black).  In  both  loading  conditions,  the  blast  wave  produces  a  rapid  rise  in  pressure  at 
around  2  ms,  which  is  qualitatively  similar  for  all  of  the  five  regions  (panels  A  and  D). 
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Figure  1.  Finite  element  results.  Pressure,  axonal  strain,  and  effective  strain  rate  measured  at  various 
locations  of  the  brain  in  both  blast  loading  scenarios. 
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What  follows  the  onset  of  the  initial  pressure  wave  are  a  number  of  smaller-amplitude  and  lower- 
frequency  pressure  waves  that  vary  from  region  to  region.  Interestingly,  the  largest  axonal  strains 
(panels  B  and  E)  occur  several  milliseconds  after  the  initial  pressure  wave.  The  large  axonal 
strains  are  the  result  of  the  angular  rotation  of  the  head.  This  motion  occurs  more  slowly  than  the 
pressure  wave  passing  through  the  head,  resulting  in  a  delay  before  the  axonal  strain  increases. 
The  effective  strain  rate  (panels  C  and  F)  shows  high-frequency  oscillations.  Strain  rate  has  been 
shown  to  contribute  to  cellular  death  in  some  brain  regions  (51)  so  it  is  important  to  consider  this 
value  along  with  the  axonal  strain.  Emphasis  should  be  placed  on  regions  that  show  a  high 
axonal  strain  and  effective  strain  rate,  such  as  the  occipital  and  temporal  lobes.  These  regions  are 
more  likely  to  show  damage,  and  it  will  be  shown  that  they  have  a  higher  network  degradation. 

The  primary  purpose  of  this  work  is  to  describe  a  new  method  for  degradation  of  the  structural 
network  and  evaluate  its  use  on  a  small  network  for  future  use  as  input  into  a  functional 
electrophysiological  model.  While  the  model  presented  in  this  work  is  currently  not  validated  for 
blast  loading  conditions,  these  simulations  provide  two  possible  scenarios  to  compare  the 
network  response.  This  allows  for  the  investigation  of  structural  network  degradation  that  results 
from  tissue  deformation  calculated  from  an  FE  simulation.  Validation  of  blast  loading  and  a 
more  detailed  description  of  the  mechanics  involved  in  the  blast  simulation  will  be  discussed  in 
future  work. 

2.3  Mapping  FE  Data  to  Voxel  Space 

To  make  use  of  neuroimaging  data,  FE  results  are  mapped  to  the  standardized  voxel  space  in 
magnetic  resonance  data,  where  each  voxel  represents  a  2-  x  2-  x  2-mm  cube  of  brain  tissue. 

This  mapping  ensures  that  the  FE  data  has  the  same  volume,  resolution,  and  Cartesian 
coordinates  as  the  diffusion  tensor  data.  This  also  allows  the  analysis  to  link  with  standardized 
brain  atlas  data  that  contains  coordinates  and  labels  for  common  brain  regions  studied  in  the 
neuroscience  field  (52).  There  are  two  possible  qualifications  where  an  element  variable  is 
assigned  to  a  voxel  as  illustrated  in  figure  2.  In  this  image,  the  black  grid  is  a  representation  of 
voxels  and  the  yellow  triangles  represent  elements.  Panel  A  shows  a  case  where  an  element  is 
assigned  to  voxels  where  the  center  of  the  voxel  is  contained  within  the  element  (red  dots).  Panel 
B  shows  an  an  element  that  does  not  contain  any  voxel  centers.  In  this  case,  the  element  is 
assigned  to  the  voxel  that  contains  the  element’s  centroid  (blue  dot).  Multiple  elements  within  a 
single  voxel  are  averaged. 
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Figure  2.  Element  to  voxel  assignment. 


2.4  Network  Construction  and  Degradation 

A  damage  parameter  based  on  cellular  death  is  then  calculated  for  each  voxel  using  an  empirical 
formula  for  cellular  death  up  to  four  days  post-insult  based  on  the  assigned  axonal  strain  and 
effective  strain  rate  (11,  50).  This  single-damage  parameter  is  then  used  as  a  means  to  degrade 
the  structural  network.  A  cell  death  threshold  for  electrophysiological  injury  is  estimated  based 
on  the  work  of  Bain  and  Meaney  (53),  who  claimed  that  electrophysiological  injury  would 
initially  occur  between  13%  and  28%  strain  in  a  guinea  pig  optic  nerve,  and  suggested  an 
optimal  threshold  of  18%  strain  at  a  strain  rate  between  30  to  60  s-i  (53).  Applying  these  values 
to  the  equation  for  cellular  death  results  in  an  average  of  2.1%  for  injury  initiation  and  3.4%  cell 
death  for  the  optimal  injury  threshold  for  this  range  of  strain  rates.  Using  these  cell  death  values 
is  beneficial  because  it  incorporates  strain,  strain  rate,  and  time  after  injury  into  a  single  value 
that  is  used  as  a  threshold  for  damage.  This  will  help  to  account  for  the  effects  of  high -rate 
loading  and  allow  for  a  prediction  of  injury  even  at  a  lower  strain  if  the  strain  rate  is  sufficient  to 
allow  the  cell  death  to  reach  threshold  values. 

The  Connectome  Mapper  Toolkit  is  used  to  parcellate  the  brain  into  83  regions  of  interest  as  in 
previous  work  (11)  based  on  the  Desikan-Killiany  brain  atlas  (52).  The  network  is  then 
resampled  to  12  cortical  ROIs  to  investigate  methods  to  maximize  compatibility  with  different 
size  networks.  These  are  the  frontal  lobe,  parietal  lobe,  occipital  lobe,  medial  temporal  lobe, 
lateral  temporal  lobe,  and  cingulate  cortex  in  both  the  left  and  right  hemisphere.  These  12  ROIs 
become  the  nodes  of  the  network.  Figure  3  shows  the  original  83-node  parcellation  from  a  lateral 
view  A  and  medial  view  B  where  colors  are  used  to  separate  the  regions.  The  smaller  cortical 
regions  in  each  lobe  are  combined  to  form  the  new  12  ROIs.  These  new  regions  are  shown  in  C 
and  D  for  the  lateral  and  medial  view,  respectively.  There  is  the  potential  for  variation  in  the 
resulting  network  measures  when  down-sampling  to  a  12-node  because  network  measures  are 
dependent  on  the  topology  of  the  network.  Networks  with  68,  1 14,  216,  446,  and  1002  nodes  are 
also  created  by  subdividing  the  original  parcellation  using  the  process  described  by  Cammoun 
et  al.  (54).  This  is  done  to  better  understand  how  this  change  in  topology  will  affect  network 
degradation. 
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Figure  3.  An  illustration  of  the  surface  parcellation.  The  original  83-node  parcellation  shown 
from  a  lateral  view  (A)  and  medial  view  (B)  is  used  to  create  a  12-node  cortical 
parcellation  shown  from  a  lateral  view  (C)  and  medial  view  (D). 

Fiber  tractography  from  diffusion  tensor  imaging  in  combination  with  the  ROIs  created  during 
the  parcellation  are  used  to  create  the  initial  network.  Fibers  that  begin  in  one  ROI  and  end  in 
another  represent  a  pathway  for  communication.  This  connection  is  represented  by  an  edge 
between  two  nodes.  Figure  4  illustrates  a  simplified  picture  where  the  grey  represents  ROIs  and 
the  black  lines  are  tractography  fibers  that  make  up  an  edge  between  the  two  ROIs.  The  squares 
of  the  grid  are  a  simplified  two-dimensional  illustration  of  voxels.  Network  edges  are  given 
weights  to  account  for  connection  strength  using  a  method  that  will  be  referred  to  as  the  fiber 
segment  method.  Each  fiber  is  divided  into  segments  based  on  the  number  of  voxels  that  it 
passes  through.  The  number  in  each  voxel  of  figure  4  is  the  number  of  fiber  segments  contained 
within.  The  weight  of  an  edge  between  two  ROIs  is  equal  to  the  total  number  of  fiber  segments 
in  all  tracts  that  connect  these  ROIs  divided  by  the  number  of  voxels  on  the  edge.  The  example 
in  figure  4A  results  in  a  weight  of  2.  This  method  is  different  from  the  weighting  method  in  Kraft 
et  al.  (11),  which  weights  edges  based  on  the  fiber  as  a  whole.  Panel  B  of  figure  4  shows  an 
example  of  a  damaged  connection.  Voxels  are  weighted  in  red  if  they  are  above  the  minimum 
cellular  death  threshold.  Numbers  shown  in  red  are  the  remaining  fiber  segments  after  damage  is 
applied.  The  weighting  of  the  edge  is  based  on  the  number  of  fiber  segments  that  remain  after 
damage  is  applied.  Dividing  by  the  number  of  voxels  in  the  edge  accounts  for  long-distance 
connections  that  would  be  composed  of  a  large  number  of  fiber  segments.  Edge  weights  are 
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degraded  by  removing  fiber  segments  from  voxels  that  meet  or  exceed  the  cell  death  threshold. 
The  percentage  of  fiber  segments  removed  is  0%  at  or  below  the  lower  threshold  of  2.1%  cell 
death  and  increases  linearly  to  100%  fiber  segment  removal  at  the  upper  threshold  of  3.4%  cell 
death. 


Figure  4.  Fiber  segment  method  edge  degradation.  An  illustration  of  an  example  edge  made  up  of  two 

regions  of  interest  connected  by  tractography  fibers.  (A)  shows  the  edge  before  degradation  and 
labels  the  number  of  fiber  segments  within  each  voxel.  (B)  shows  an  example  of  cell  death  data 
applied  to  degrade  the  edge. 

Weights  are  resampled  to  a  Gaussian  distribution  N(.5.\  )  into  the  range  [0,1].  This  is  based  on 
procedures  developed  by  Honey  et  al.  (55),  who  stated  the  assumption  that  interregional 
physiological  efficacies  would  not  span  the  large  range  that  is  seen  in  the  data  before  resampling. 
This  Gaussian  weight  degradation  is  defined  as  G  =  Go[WAVo],  where  G  is  the  damaged 
Gaussian  weight,  Go  is  the  undamaged  Gaussian  weight,  W  is  the  damaged  weight,  and  Wo  is  the 
undamaged  weight. 

The  fiber  segment  method  is  beneficial  because  it  avoids  the  complete  removal  of  edges  seen  in 
previous  work  by  Kraft  et  al.  (11).  This  method  results  in  a  lower  estimation  of  network  loss 
compared  to  previous  methods  where  the  entire  fiber  is  removed  but  the  pattern  of  degradation 
remains  the  same.  However,  there  are  several  potential  limitations  of  this  weighting  method. 

First  is  the  assumption  that  fibers  that  pass  through  a  region  with  a  high  damage  value  retain 
some  capacity  for  communications.  While  this  would  not  be  valid  for  a  single  axon,  it  is 
reasonable  in  the  case  of  a  fiber  segment  representing  the  volume  of  a  voxel  at  a  resolution  of 
2  mm3.  This  segment  of  fiber  will  be  composed  of  a  large  number  of  axons,  some  of  which  may 
be  damaged  while  others  survive  and  can  propagate  an  electrophysiological  signal.  Next,  due  to 
limited  experimental  data,  this  method  cannot  represent  smaller-scale  injury  mechanisms  and 
cannot  provide  a  realistic  quantitative  value  of  axonal  death.  The  use  of  a  threshold  limits  the 
range  of  damage  prediction  to  this  value.  This  is  useful  because  it  allows  for  greater  sensitivity 
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for  detecting  minor  injury,  but  reductions  in  edge  weight  do  not  necessarily  correspond  to  the 
same  reduction  in  biological  signaling  capacity  of  cells  that  connect  the  ROIs.  Despite  these 
issues,  this  method  is  useful  because  it  offers  an  estimate  of  the  scale  of  damage.  It  is  used  to 
determine  where  damage  occurs  and  the  relative  severity  of  the  damage  seen  in  different  regions 
of  the  brain.  This  will  be  helpful  for  differentiating  possible  injury  based  on  different  loading 
conditions.  This  model  will  be  improved  to  produce  a  more  biologically  accurate  prediction  of 
structural  network  damage  as  relevant  experimental  data  becomes  available. 


3.  Results  and  Discussion 


3.1  Network  Degradation 

Figure  5  shows  the  fiber  tractography  for  frontal  loading  in  panels  A-D  and  side  loading  in 
panels  E-H  for  the  24-,  48-,  72-,  and  96-h  time  points.  The  anterior  direction  is  to  the  right  and 
the  posterior  direction  is  to  the  left.  Undamaged  fibers  are  colored  based  on  their  direction,  with 
green,  blue,  and  red  representing  anterior-posterior,  inferior- superior,  and  left-right,  respectively. 
Fiber  segments  that  are  removed  from  the  model  are  displayed  in  black.  Relatively  few  fibers  are 
removed  in  the  frontal-loading  condition  until  72  h.  At  96  h,  the  location  of  removed  fibers  is 
distributed  throughout  the  volume  of  the  brain.  The  side  loading  shows  a  high  concentration  of 
fiber  segments  that  are  removed  in  the  occipital  lobe  starting  at  48  h,  which  corresponds  to  the 
high  strain  and  strain  rate  that  is  seen  in  the  FE  time-history  plots  in  figure  6.  The  location  of 
removed  fibers  continues  to  expand  up  to  96  h.  This  image  format,  developed  by  Alper  et  al. 

(56),  displays  the  undamaged  edge  weight  colored  by  edge  strength  overlapped  with  a  smaller 
square  that  is  colored  by  the  damaged  edge  weight.  Because  this  is  an  undirected  network  that 
produces  a  symmetrical  matrix,  the  two  halves  are  weighted  to  show  different  values.  The  blue- 
weighted  half  shows  the  edge  strength  and  the  red- weighted  half  shows  the  normalized  edge 
strength.  The  regions  of  interest  are  the  medial  temporal,  lateral  temporal,  frontal,  parietal, 
occipital,  and  cingulate  cortices  in  the  right  hemisphere  for  regions  1-6  and  left  hemisphere  for 
regions  7-12. 
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Figure  5.  Evolution  of  fiber  segment  removal,  using  empirically  based  cellular  death  predictions 

obtained  from  in  vitro  models  of  neural  tissues.  Local  strain  and  strain-rate  values  computed 
from  FE  simulations  are  used  to  specify  injury. 

Basic  measurements  are  calculated  from  the  data  represented  in  figure  6.  The  weight  at  each 
edge  of  the  network  is  represented  by  a  single  element  of  the  matrix  that  corresponds  to  a 
connection  between  two  specific  ROIs.  This  is  the  most  localized  measure  to  investigate  damage 
between  two  regions.  The  degree  is  calculated  as  the  sum  of  matrix  elements  along  a  single 
column  or  row  for  any  one  node.  This  nodal  measure  shows  damage  to  any  one  node  based  on  all 
of  its  connections.  The  inner  square  of  each  matrix  element  represents  the  value  at  the  current 
time  point  for  comparison  to  the  outer  square  that  represents  the  initial  value.  The  regions  of 
interest  are  medial  temporal,  lateral  temporal,  frontal,  parietal,  occipital,  and  cingulate  cortices  in 
the  right  hemisphere  for  regions  1-6  and  left  hemisphere  for  regions  7-12.  The  red-weighted  half 
of  the  matrix  shows  the  normalized  edge  strength,  and  the  sum  of  all  edge  weights  is  calculated 
by  the  sum  of  all  matrix  elements  in  the  blue- weighted  half  of  the  matrix.  This  value  provides  a 
global  account  of  the  network  to  show  how  it  has  changed  as  a  whole. 
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Figure  6.  Undirected  connectivity  matrices. 

The  frontal-loading  condition  is  shown  in  figure  6,  panel  A.  This  figure  shows  that  there  is  little 
damage  applied  to  the  network  on  a  global  level  with  1.8%  reduction  in  the  total  edge  strength  at 
96  h.  Because  of  the  subtle  change  in  this  global  measure,  it  is  important  to  consider  localized 
measures  of  damage.  Table  1  shows  the  percentage  change  in  both  degree  and  edge  weight  of  the 
three  most  damaged  regions  or  connections  in  either  case. 

Table  1.  Percent  change  in  degree  and  edge  weight  for  the  most  damaged  regions  in  the  frontal- 
loading  condition. 


Region 

(No.) 

Degree 

(%) 

Edge  Connection  Number 
(No.) 

Weight 

Reduction 

(%) 

Cingulate  R  (6) 

3.4 

Parietal  R  (4) 

Med.  Temporal  L  (7) 

11.3 

Parietal  R  (4) 

2.6 

Lat.  Temporal  R  (2) 

Frontal  L  (9) 

5.7 

Cingulate  L  (12) 

2.3 

Parietal  R  (4) 

Cingulate  R  (6) 

5.5 

The  side-loading  condition  is  shown  in  figure  6  panel  B.  A  larger  global  damage  is  seen  in  this 
case  where  the  percent  change  in  total  edge  strength  is  10.3%  at  96  h.  Many  connections  show 
damage  resulting  in  larger  change  in  the  degree  and  edge  weight.  These  two  measures  are  shown 
as  a  percentage  change  in  table  2  for  the  most  damaged  regions  or  connections.  While  these 
results  reconfirm  that  regions  with  high  strain  and  strain  rate  seen  in  the  FE  results  show  higher 
network  degradation,  they  also  show  damage  to  some  additional  ROIs.  The  right  parietal  lobe 
had  low  strain  compared  to  other  regions  in  the  orange  trace  of  figure  1,  panel  E. 
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However,  the  right  parietal  region  is  connected  to  other  regions  with  high  damage.  As  a  result, 
the  parietal  lobe  is  at  risk  to  be  affected  by  damage  to  white-matter  pathways  from  these  distant 
brain  regions. 

Table  2.  Percent  change  in  degree  and  edge  weight  for  the  most  damaged  regions  in  the 
side-loading  condition. 


Region 

(No.) 

Degree 

(%) 

Edge  Connection 
(No.) 

Weight 

Reduction 

(%) 

Occipital  R  (5) 

20.1 

Lat.  Temporal  R  (2) 

Parietal  R  (4) 

35.5 

Lat.  Temporal  R  (2) 

16.4 

Occipital  R  (5) 

Parietal  L  (10) 

28.4 

Parietal  R  (4) 

15.5 

Parietal  R  (4) 

Occipital  L  (1 1) 

28.3 

The  betweenness  indicates  the  presence  of  a  hub  in  the  network.  This  measure  is  calculated  for 
the  initial  network  to  gain  a  better  understanding  of  which  ROI  would  potentially  create  a  more 
widespread  effect  on  the  network  if  damaged.  Figure  7  shows  the  betweenness  for  all  regions. 
The  highest  values  are  seen  in  the  parietal  region  in  both  hemispheres  followed  by  the  frontal 
lobe  of  both  hemispheres.  The  right  parietal  region  is  among  the  top  three  most  affected  regions 
for  degree  in  tables  1  and  2.  This  indicates  that  short  paths  traversing  the  parietal  region  may  be 
especially  disrupted  as  a  result  of  damage  to  the  structural  network. 


Figure  7.  The  betweenness  value  for  each  ROI.  This  plot  shows  that  the  parietal 
and  frontal  regions  allow  for  the  greatest  number  of  multiedge 
connections  along  paths  that  minimize  edge  cost. 
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3.2  Network  Measure  Scale  Invariance 

An  assessment  of  network  measures  is  performed  on  68-,  114-,  216-,  446-,  and  1002-node 
cortical  networks  in  addition  to  the  12-node  network  to  better  understand  how  the  number  of 
ROIs  affects  how  damage  is  applied  to  the  network.  These  additional  networks  are  the  result  of 
subdividing  the  original  parcellation  (54).  There  is  great  variation  in  the  value  of  network 
measures  computed  on  the  initial  undamaged  network  at  these  different  parcellation  levels.  This 
is  shown  for  a  selective  representation  of  measures  in  figure  8,  panels  A-D.  In  panel  A,  the  total 
edge  strength  of  the  network  increases  with  the  number  of  nodes  because  the  networks  composed 
of  a  larger  number  of  nodes  will  have  more  edges.  Both  global  efficiency  and  clustering 
coefficients  (B  and  C)  show  significant  change  between  12  and  114  nodes.  There  is  less 
variability  in  these  values  between  1 14  and  1002  nodes.  The  initial  high  values  are  the  result  of 
the  12-node  network  being  the  closest  to  a  complete  network  where  all  regions  are  connected. 
The  network  topology  is  similar  to  a  single  cluster  in  the  small  networks  rather  than  the 
collection  of  localized  clusters  that  is  seen  in  the  larger  networks.  The  efficiency  and  clustering 
coefficient  are  reduced  as  the  network  connections  become  more  sparse  after  additional  nodes 
are  used  to  form  the  network.  The  increase  in  modularity  in  panel  D  at  high  parcellations 
indicates  the  formation  of  localized  highly  connected  regions  that  are  not  present  in  the  small 
networks. 

It  is  important  to  consider  normalized  network  measures  due  to  the  large  variability  of  measures 
applied  to  networks  of  different  scales.  The  normalized  network  measure  is  defined  as 
N  =  [M/Mo]  100%,  where  N  is  the  normalized  measure,  M  is  the  damaged  measure  and  Mo  is 
the  undamaged  measure.  By  calculating  this  normalized  measure,  the  amount  of  damage  is  more 
consistent  across  all  parcellations.  Panel  E  of  figure  8  shows  the  normalized  measure  of  total 
edge  strength,  global  efficiency,  clustering  coefficient,  average  shortest  path,  assortativity,  and 
modularity  at  96  h.  All  networks  shows  similar  normalized  values  of  total  edge  strength,  global 
efficiency,  clustering  coefficient,  and  average  shortest  path.  Of  these  measures,  the  total  edge 
strength  in  the  side-loading  condition  has  the  greatest  variation  with  at  most  2.05%  difference 
between  networks  with  any  number  of  nodes.  However,  modularity  and  assortativity  show 
greater  deviation  at  the  12-node  network  compared  with  other  parcellations. 
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Figure  8.  Effects  of  network  scale. 


3.3  Network  Measure  Sensitivity 

In  addition  to  being  scale-invariant,  a  network  measure  should  be  sensitive  to  different  loading 
conditions  and  be  able  to  predict  different  magnitudes  and  locations  of  damage.  Efficiency, 
clustering  coefficient,  and  average  shortest  path  maintain  a  consistent  percentage  of  damage 
across  all  parcellations.  These  measures  are  chosen  for  a  more  detailed  analysis  of  sensitivity. 
The  efficiency  is  useful,  as  it  can  be  calculated  on  a  global  and  local  scale.  The  ratio  of  average 
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clustering  coefficient  and  average  shortest  path  provides  a  measure  that  is  comparable  to  the 
global  efficiency  (35).  Figure  9  shows  the  relative  total  edge  weight,  efficiency,  and  ratio  of 
average  clustering  coefficient  and  average  shortest  path  through  96  h.  The  normalized  values  of 
these  measures  are  98.1%,  98.0%,  and  96.1%,  respectively,  for  the  frontal  loading  and  89.7%, 
89.6%,  and  88.0%,  respectively,  for  the  side  loading.  All  three  measures  are  reduced  consistently 
from  0  to  96  h.  For  the  12-node  network  and  this  level  of  damage,  the  more  complex  network 
measurements  such  as  efficiency  do  not  provide  additional  information. 


Figure  9.  Normalized  network  measures.  For  this  12-node  network,  all  measures 
show  a  similar  trend.  Note  that  side  loading  resulted  in  a  larger  change 
in  these  values. 

An  analysis  of  network  measures  is  performed  to  compare  the  strength  of  connections  between 
ROIs  within  the  same  hemisphere  with  connections  that  cross  hemispheres  to  further 
differentiate  between  loading  conditions.  To  do  this,  nodes  are  removed  from  one  hemisphere  to 
calculate  the  global  efficiency  of  the  remaining  nodes.  This  results  in  a  measure  of  global 
efficiency  for  the  right  hemisphere  only  and  for  the  left  hemisphere  only.  Figure  10  shows  the 
normalized  global  efficiency  for  the  right  hemisphere,  left  hemisphere,  and  the  whole  brain  for 
times  up  to  96  h.  This  shows  that  the  normalized  efficiency  has  a  9.0%  difference  in  the  right 
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hemisphere  compared  with  the  left  hemisphere  for  the  side-loading  condition  at  96  h.  Note  that 
the  global  efficiency  of  the  brain  as  a  whole  in  the  right  side-loading  condition  is  reduced  at  a 
level  similar  to  the  right  hemisphere.  This  indicates  that  the  damage  to  edges  that  connect  regions 
across  hemispheres  is  great  enough  to  cause  more  widespread  effects.  The  frontal  loading 
produces  a  more  consistent  reduction  in  normalized  global  efficiency  for  both  hemispheres  with 
only  a  1.3%  difference  between  the  left  and  right.  These  results  are  expected  based  on  the 
location  of  fiber  segments  removed.  This  shows  that  network  measures  are  useful  to  quantify  the 
difference  in  damage  to  the  two  hemispheres  of  the  brain  based  on  the  loading. 
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Figure  10.  Lateral  comparison  of  efficiency. 

3.4  Discussion  and  Future  Work 

There  are  two  key  goals  described  in  this  research.  The  first  is  to  evaluate  the  possibility  of  using 
a  network  with  a  low  number  of  nodes  by  investigating  which  network  measures  are  scale- 
invariant.  Next  is  an  assessment  of  the  sensitivity  of  network  measures  for  different  loading 
conditions.  The  achievement  of  these  two  goals  provides  a  better  understanding  of  how  the 
newly  developed  weighting  methods  produce  damage  in  the  structural  network  of  the  brain. 
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Results  suggest  that  some  graph  measures  are  not  scale-invariant.  Topological  information  is 
distorted  when  down-sampling  to  a  12-node  network.  As  a  result,  measures  such  as  assortativity 
and  modularity  are  not  useful  for  this  small-scale  network.  However,  the  normalized  measure  of 
damage  is  more  consistent  for  the  global  efficiency,  total  edge  weight,  clustering  coefficient,  and 
average  shortest  path  across  networks  with  12  to  1002  nodes.  This  shows  that  these  network 
measures  will  be  useful  for  comparing  a  network  of  any  size. 

Our  results  indicate  that  the  new  weighting  method  produces  different  patterns  of  damage  to  the 
structural  network  corresponding  to  the  different  loading  conditions  presented  here.  Network 
measures  such  as  global  efficiency,  total  edge  weight,  clustering  coefficient,  and  average  shortest 
path  are  used  to  show  the  difference  in  damage  applied  to  the  network  on  a  global  level. 
However,  all  of  these  measures  predict  a  similar  magnitude  of  global  damage.  On  a  12-node 
network  for  these  loading  conditions,  there  is  no  advantage  in  using  more  complex  network 
measures.  This  shows  that  local  measures  are  necessary  to  understand  the  how  damage  is  applied 
to  the  network.  Nodes  on  the  side  closest  to  the  blast  origin  were  more  damaged  in  the  side¬ 
loading  condition.  However,  there  is  only  a  small  amount  of  variation  in  the  degree  in  the  front- 
loading  condition.  For  this  loading  condition,  individual  edge  weights  are  the  most  effective 
indication  of  where  damage  takes  place.  While  this  is  a  very  simple  measure,  it  is  the  most 
reliable  estimation  of  damage  for  a  small  network  such  as  the  12-node  network  investigated  in 
this  research.  In  addition,  individual  edge  weights  show  the  greatest  sensitivity  to  minor  damage 
when  relatively  few  fiber  segments  are  removed  from  the  model.  Because  of  this  sensitivity, 
changes  in  individual  edge  weights  are  suggested  for  use  to  determine  changes  in  coupling 
strength  in  an  electrophysiological  model. 

In  future  work,  the  ability  to  make  use  of  a  functional  brain  network  model  based  on  a  degraded 
structural  network  will  allow  for  additional  insight  into  the  possible  effects  of  blast  loading  to  the 
brain.  In  the  electrophysiological  model,  the  simulated  brain  region  nodes  oscillate  at  particular 
intrinsic  frequencies,  and  functional  networks  are  formed  by  connections,  or  edges,  between  the 
nodes.  The  edge  weighting  in  an  electrophysiological  model  will  be  determined  by  the  edge 
weights  estimated  in  multiple  degraded  structural  networks  at  different  time  points  (e.g., 
undamaged  and  at  24,  48,  72,  and  96  h).  This  will  allow  for  the  adjustment  of  coupling  strength 
between  simulated  brain  nodes  to  represent  the  functional  aspect  of  the  network,  revealing  how 
changes  in  edge  weighting  based  on  the  loading  condition  manifests  in  changes  in  the  functional 
brain  network.  Investigation  of  instability  in  the  functional  network  due  to  simulated  structural 
network  damage  predicated  from  an  FE  simulation  can  help  improve  not  only  our  understanding 
of  the  physical  damage  to  the  brain  but  also  the  possibility  of  electrophysiological  impairment. 
This  provides  a  technical  framework  to  use  computational  methods  to  understand  how  structural 
injury  links  to  abnormal  electrical  oscillations  that  disrupt  the  functional  connectivity  patterns 
associated  with  successful  task  performance.  In  short,  this  framework  links  structural  damage 
with  disrupted  functional  connectivity  that  may  underlie  cognitive  impairment. 
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4.  Conclusions 


Our  larger  research  program  on  brain  structure-function  couplings  at  ARL  integrates  modeling 
and  experimentation  to  elucidate  the  links  among  structural  connectivity,  functional  connectivity, 
and  behavioral  performance.  The  research  described  in  this  report  details  a  technical  approach  to 
improve  the  understanding  of  neurotrauma  through  the  combined  use  of  the  finite  element 
method,  structural  neuroimaging  methods,  and  graph  theory  approaches  from  network  science. 
This  report  presents  an  expansion  of  our  previous  methods  using  a  physics-based  simulation  to 
inform  structural  brain  network  degradation.  Our  new  weighting  method  using  fiber  tract 
segments  provides  increased  sensitivity  to  structural  network  damage  based  on  graded  levels  of 
cellular  tissue  damage.  Our  results  indicate  that  the  new  weighting  method  can  differentiate 
between  loading  conditions  in  the  simulated  damage  to  the  structural  network  model.  The  results 
also  indicate  that  some  graph  measures  are  not  scale-invariant  and  should  be  selected  based  on 
the  size  of  the  network.  Evaluation  of  a  12-node  network  shows  that  the  most  reliable  measure  of 
degradation  is  a  reduction  in  individual  edge  weights. 

Future  work  will  use  time-evolving  estimates  of  degraded  edge  weights  to  inform  the 
connections  between  simulated  brain  nodes  in  an  electrophysiological  model  and  examine  how 
different  blast-loading  conditions  change  electrical  activity  in  a  functional  brain  network.  By 
understanding  which  brain  regions  see  diminished  capacity  for  communication  in  the  structural 
network,  it  may  be  possible  understand  instabilities  in  functional  connectivity  that  may  underlie 
the  cognitive  deficits  characteristic  of  traumatic  brain  injuries.  Substantial  experimental  work  is 
needed  to  better  understand  how  to  relate  the  network  degradation  described  in  this  work  to 
biological  injury.  However,  this  method  provides  a  preliminary  pathway  to  construct  more 
realistic  models  of  mechanical  tissue  deformation  using  information  about  white-matter  tracts 
from  diffusion  tensor  imaging  and  to  develop  tools  to  understand  how  structural  injury  translates 
to  functional  network  changes. 
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